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Abstract 

Adjoint field methods are both elegant and efficient for calculating sensitivity infor- 
mation required across a wide range of physics-based inverse problems. Here we provide 
a unified approach to the derivation of such methods for problems whose physics are pro- 
vided by Poisson's equation. Unlike existing approaches in the literature, we consider in 
detail and explicitly the role of general boundary conditions in the derivation of the asso- 
ciated adjoint field-based sensitivities. We highlight the relationship between the adjoint 
field computations required for both gradient decent and Gauss-Newton approaches to 
image formation. Our derivation is based on standard results from vector calculus cou- 
pled with transparent manipulation of the underlying partial different equations thereby 
making the concepts employed here easily adaptable to other systems of interest. 

1 Introduction 

Geophysical imaging modalities based on the inversion of the Poisson's equation include elec- 
trical resistance tomography (ERT) [l], electrical impedance tomography (EIT) and in- 
duced polarization (IP) |3|. Additionally, the same physical model underlies electrical capac- 
itive tomography (ECT), employed for nondestructive evaluation f4','5] as well as aspects of 
the diffuse optical tomography inverse problem arising in brain and breast imaging f6^. Most 
all inversion methods require the calculation of sensitivity information as part of the imaging 
algorithms; i.e. the functional derivative of the cost function or a part of the cost function 
with respect to the unknown physical quantity being imaged. 

Adjoint field methods represent an analytically elegant as well as computationally at- 
tractive approach for sensitivity calculations and have been considered for various imaging 
problems. Among the main contributions in this area we highlight the work by Somersalo 
et al. in which the adjoint calculations are analyzed for the EIT problem using weak forms 
of the Poisson's equation [T]. In the context of electromagnetic imaging, Dorn et al. in [s] 
have considered reconstruction of the complex conductivity using magnetic measurements, for 
which they use the Maxwell's equations and derive the required sensitivity information. In 
an alternative contribution, shape based reconstruction of the imaging domain is considered 
using the differential form of the Helmholtz equation as the governing modality [o]. Inverse 
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problems based on the integral forms of the Helmholtz equation and full wave models taking 
into account the incident and scattered field components are another challenging problem for 
which adjoint methods play a significant role in simplifying the sensitivity calculations. For 
more details about the application of adjoint technique in these problems, an interested reader 
is referred to some of the work of Abubakar et al. mainly in the context of contrast source 
inversion (CSI) fTo}{l2 . 



While the adjoint field approach has been applied to Poisson inverse problems, the details 
of its use depend quite heavily on the algorithmic method being employed by the inverse 



solver. More specifically, the use of a gradient decent-type of approach [13-15 requires a 
different adjoint calculation than a Gauss- Newton method for which a full Jacobian must be 
determined [16| - |18| . Additionally, it turns out that given the latter, the former can easily be 
obtained. 

It is certainly true that the literature contains a number of similar derivations based on 
operator-theoretic principals 14 ,19 , employing weak formulations [7 20 



21 
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, or using physical 



24 . While the final, 



concepts such as power conservation [17] and the reciprocity theorem 
analytical expressions for the sensitivity information are quite similar (if not the same) to that 
which is derived here, the analytical methods used in the various derivations differ markedly 
in terms of brevity, clarity, and the level of mathematical background underlying the analysis. 
Moreover, as noted previously, the adjoint computations can be put to a number of different 
uses depending on the nature of the optimization technique being employed to solve the 
inverse problem; i.e., gradient decent or Newton- type. The differences in these uses are not 
readily apparent in the adjoint-field derivations for the Poisson problem currently available 
in the literature. 

The primary contribution of this paper is the presentation of a clear and detailed derivation 
based on easily accessible, vector calculus identities (following one of Norton's methods in 
|25| where the Helmholtz equation was considered) of adjoint based sensitivity calculations 
for Poisson based inverse problems. We highlight the mathematical relationship between 
the gradient descent and Newton- type adjoint forms and consider problems with general 
type of boundary conditions, including Neumann, Dirichlet or mixed boundary conditions. 
While we certainly acknowledge that some of the results presented in the paper are well 
known, we feel that the pedagogical approach we have taken, our comprehensive and detailed 
treatment of boundary conditions, and the explicit treatment of the calculations required 
for gradient decent and Newton-type uses of the approach represent a contribution to the 
literature. Moreover, by employing a minimum of mathematical abstraction and a step- by- 
step derivation of all results, we hope that the approach we take here can provide the interested 
reader with the tools needed to apply the ideas to problems governed by other physical models. 



2 General Problem Formulation 

For an imaging domain with surface F, the forward model of interest motivated by several 
geophysical applications is 

V-(o-V0)=s in J], (1) 

ao-|^ + ^0 = O onF. (2) 
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Figure 1: An ERT application: reconstruction of the subsurface conductivity distribution 
based on the potential measurements at inserted boreholes. A Neumann boundary condition 
(B.C.) considered for the top surface and a mixed or Dirichlet type B.C. considered for the 
remaining boundaries 



In the context of electrical tomography as a basic modality, (f) is the electrical potential, s 
current source and a the conductivity. All quantities are assumed to be functions of three 
spatial coordinates, r = [x y z]^. Regarding the boundary condition ([2|, the notation d(f)/dii 
denotes the normal component of the gradient of (j) on F, and a(r) and /3(r) are functions 
defined on the surface F (i.e., r e F) which are not simultaneously zero. To maintain generality 
and consider problems with different types of boundary conditions on different regions of F, 
we do not make any continuity assumptions about a and f3. Also, generalization of the results 
to nonhomogeneous boundary condition will be accomplished later in this letter. 

A common geophysical problem associated with the model in ([T])-([2]) is shown in Fig. 1. 
In this problem F consist of F„, the interface between the earth and air where we impose 
a zero current condition (/3 = 0) and Tm, a surface where the values a and /3 are chosen to 
model an infinite half-space |26| . As an approximation to an infinite half-space, the boundary 
condition on F^ may be replaced with a zero potential condition (a = 0) when F^ is far from 
the sources of current f27j. Reconstruction of a based on the measurements of (f) at some 
points in the domain is the goal of this tomography problem. 

For simplicity here, we consider only real-valued conductivity (that is the ERT problems) 
although the approach can easily be adapted for use in estimation of complex valued con- 
ductivity as is encountered in EIT as well as conductivity/chargability as are desired in IP 
experiments. To keep the notation simple, we also consider the case where data are collected 
from a single source of current. More generally, in a tomographic imaging problem one would 
illuminate with many sources Sp for p - 1, 2, P. 

Central to the derivation that follows is the impact of a small change in conductivity to 
the system ([Tll-Q. Consider a perturbation to the conductivity, a a + 6a resulting in 
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(j)^ (f) + 5(f). Using these in ([Tf-Q gives 

V- ((cr + (5cr)V((/> + 50)) = s onO, (3) 

a{a + 5a)-^{(t) + 54>) + ^{(t) + 6(t))^Q on T. (4) 

an 

Now, expanding Q as 

V • (crV(/)) + V • {5aV4>) + V • (crV(50) + V • {SaVScj)) = s, 
keeping terms of first order and using ([l]) in ([3]) and ([2]) in (|4]) yields 

V ■ {6aV<t>) + V ■ {aV6(j)) ^0, on 0, (5) 

a{a-^d(l) + 5a^) + f56())^0 onT. (6) 

an an 

The goal of the inverse problem is to estimate a from observations of the potential collected 
at a finite set of points in space rj, i - 1,2,---,N. The estimate is obtained by minimizing the 
least squares cost function 

1 ^ 2 

^ j=l 

= ^-B{afBia), (7) 

where 4>obsi^i) are the observations of potential at location rj and (j){ri) is implicitly a function 
of a through Q-Q. The length column vector E(cj) contains the residuals, ei{a) = 
(/>(rj) - 4>obs{T^i) for i - 1,2,---,N. A few remarks are in order regarding ([T]): 

• The use of a least squares formulation is not terribly restrictive. For example, a one- 
norm type of cost function could be accommodated by populating E with a suitably 
smoothed version of \4>{vi) - (/'o6s(ri)l^^^- Though somewhat tedious, the results in this 
letter could then be generalized. 

• Generalization of ^ to account for weighting of the residuals is also not difficult. To 
avoid the resulting notational burdens, we choose to not include these details. 

• In many cases, the cost function includes both a data misfit term as in ^ as well as 
additional regularization. We refer the reader to [28] for details concerning the treatment 
of Tikhonov-type regularization schemes. The adjoint field calculations of interest here 
pertain only to the data term. 

As discussed in Section[l| two type of sensitivities are of interest depending on the nature of 
the inversion algorithm. Methods based on gradient information such as the steepest decent, 
require the functional derivative of the cost function with respect to a. More formally, a 
linear integral operator is sought which maps small perturbations of the conductivity, 5a, to 
corresponding changes in J. In the case of Newton type methods and more specifically the 
Gauss-Newton approach, one requires A^ linear operators relating 5a to perturbations in the 
individual residuals, 5ei. We start with the latter. Subsequently we indicate the changes 
to the derivation required to obtain the former and also how the adjoint calculations used 
in the Newton type approaches can be assembled to obtain the gradient-based sensitivity 
information. 
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3 Adjoint Field Calculations 

To start, we consider the variation in a due to a small change in a. We have 

5ei{a) = J0(r)|r=r,, 
which can be written volume integral 

5ei{a) = j^5{r-Vi)5(t){v)dr, (8) 

where 5{.) denotes the Dirac delta function. The objective here is to find the linear integral 
operator relating 5ei to 5a. Toward this end we define the i-th adjoint source Si as 

Si{r) = 5{r-Vi). 

It is now useful to define the potential, as the solution to the adjoint system 

V • (o-V(^i) = Si on il, (9) 

aa?^ +/30i = onT. (10) 

an 

Physically, (fn is the potential field arising from the solution to the adjoint Poisson's equation 
(which in this case is the same as the original since Poisson's equation is self-adjoint) due to 
a point source of amplitude unity located at the i-th receiver. From ^ and ^ we conclude 
that the perturbation to the residuals can be written in terms of the adjoint field as 

Sei{a) ^ J^V ■{a\/4>i)S4>{r)dr. (11) 

The goal is to find a collection of kernel functions, ki(r) such that (11) can be expressed as 

5ei{a)^ j^ki{r)5a{v)dr. (12) 



The resulting kernels are known as the Frechet derivative of Cj with respect to a [29] . 

To find these kernels, we make extensive use of the following identity derived from Green's 
theorem [30| for vector function F and scalar function g 

J^F-Vgdr+ J^g{V ■F)dr^ J^gF-dS, (13) 

where / \^ ■ dS - f ( V • n) dS denotes the surface integral over F of the normal component of 
the vector field V. 

We begin by taking g - 5(j) and F = crV^i in ( 11 ) to obtain 

5ei(cr) = - f a(S/ci>i)-(S/6(j))dr+ f a\/4)i6(j) ■ dS. (14) 
Jn Jr 

Next using g - (pi and F = aV5(/) in the first term on the right hand side of (14), we have 

Sei{a) = (j)iV ■ ((tV50) dr 

- J^(y{VH)<Pi-dS + J^avij)i6(i)-dS. (15) 
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Prom ([5]), V • 6a\/(f) = -V • aS/5(f) which we use in the first term on the right hand side of (15) 
to arrive at 



I dr 



5ei{a) = - ^ 0iV • {SaV(f>) ^ 

- a{\/6(p)(t)i -(18 + a{V(pi)d(j) ■ dS. 



Appeahng once more to (13) with g - (pi and F = 5a\I(j) in the first term of (16) gives 



(16) 



(17) 



We now rewrite the surface integral term on the right hand side of ( |17[ ) as the integration 
over the normal component 



5ei{(T) = 6a{S/(l)i) ■ (V(/>) dr 

- (^a{S/6(t>)(f>i + 6a{S/(l))4>i - cr( V^j)'^'/') ' dS. 



-f(4 

Jr V an 



+ 5a—(j)i-a 
on a 



^5(p)dS, 
n / 



(18) 



and show that it is zero. For this purpose we multiply both sides of (10) by to arrive at 



P6(l)(f>i + aa—^5(j) - 0. 
an 

Using ([6]) to replace the term in ( 19 ) results in 



/ d 

\ an 



+ 5cr— 0i - a—^6(t)) = 0, 
an an ' 



on r. 



(19) 



(20) 



The expression within the brackets in (20) is the same as the integrand in (18). Considering 



this term, if a it Q for F^ c F, then clearly the inside bracket expression becomes zero on F^. 
For the remaining surface F \ Fq, that a = 0, we certainly have /3 ^ since a and /? may not 
be simultaneously zero and using this fact in m\ and ( 10 ) would result in 50 = and 0j = 



which again make the inside bracket term zero. Therefore for all r e F the inside bracket term 



is zero and hence the surface integral in (18) vanishes, resulting in 



5ei{a) = 6a{S/(t>i) ■ ( V0) dr. 



(21) 



This result expresses the perturbation to as a 



accordingly identifying the kernel function in (16) as ki{r) - (V(^i(r)) • (v0(r)). 
In case of a system with nonhomogeneous boundary condition 



inear operator acting on 6a as desired, and 



V • ((tV0) = s 
on 



in $7, 
on F, 



(22) 
(23) 



perturbing cr a+6a and 0-5- (l)+6(j) would result in the same equations as ([5])-([6| and ^ will be 
cancelled. Based on this result, by solving the same forward problem as ([9|-(|10[) and finding 
(pi, identical forms of sensitivity as (21) will be obtained for the case of nonhomogeneous 
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boundary condition. However, we should note that although the sensitivity results are valid, 
by definition of an adjoint system, ([9])-(10) is not the adjoint of the nonhomogeneous system 
(22)~(23). More details in this regard are provided in the Appendix. 

In cases where gradient decent-type optimization methods are used for imaging, one re- 
quires the functional derivative of the cost function, J{cr) in ([T]) with respect to the conduc- 
tivity. The required variation now is 

N 

i=l 

which similarly may be rewritten in an integral form as 

N ^ 
SJ{'T)-Zi^(^^)-^obs{ri)) / 6{r-r,)6<P{r)dr 

i=i -'^ 

Ti) - 0ofe(ri))(5(r - r,;))(50(r) dr. 



N 



There are two ways to determine the resulting Frechet derivative. On the one hand, if we 
define the composite adjoint source 

N 

sir) = {<Pi^^) - 4>obsiri))6ir - r,) (24) 

i=l 

the derivation provided above follows unaltered and we conclude 

SJ{a)^ J^Sa{V4>)-{V^)dr, (25) 



where now (pis a. single adjoint field computed according to ([9j)-( 10 ) with s replacing s,. Thus, 
for gradient decent methods, one requires only two solutions of Poisson's equation: one for 
the source and one for the composite adjoint source as opposed to 1 + solves needed for 
a Gauss-Newton scheme. Alternatively, since s - J^i {(^i^i) ~ 4>obs{T^i))si, by the linearity of 
Poisson's equation and (24) 

N 

i=l 
N 



= ^((/)(ri)-0ofe.(ri)) f 5a{V^i)-{V(t))dr, 



i=l 

which is the desired link between the two uses for the adjoint field in sensitivity calculations. 

As a concrete example, consider a 3D imaging problem where cr(r) is discretized to form 
a vector cr = [cji, ctm ]"^- Here the j-th element corresponds to a voxel Qj in the domain, 
over which the conductivity is assumed to be uniformly aj. To iteratively minimize the least 
squares problem ^ using a gradient descent approach, the A:-th iteration to find cr^^'^^') is 
performed by moving along the negative direction of V J as 

^(fc+i) _^(fc) ^ _^ik)^j{k) 

._^W(ve(^'))^e(^) 
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where 7*-*^^ > is a step size 31 . On the other hand, for Newton type methods, the updating 
is performed through solving the fohowing system for ct*-*^'*'^^ 



g(fe)^^(fc+i) 



where B is a positive definite approximation to the Hessian of J, such as (VE)"^VE in a 
Guass-Newton approach or (VE)^VE + 7I in a Levenberg-Marquardt case jsij. Clearly, in 
using a gradient descent approach only VJ (basically the product of VE-^E) is required at 
every iteration while in a Gauss-Newton approach VE itself is required. The (i, j)-th element 



of the N X M matrix VE is dei/daj. By taking 6a in (21) to be dajXUji^) with XA the 
indicator function over a set A, we have 

Bp r 

p= / (V</..)-(V0)dr. 

Thus, the cost to calculate VE is solving one forward problem to compute (j), N adjoint 
problems to determine (pi, and, for each i the evaluation of M integrals to determine each row 
of the matrix. However, in a gradient descent approach where only the vector V J of length 



is required each element may be calculated through (25) as 



dJ 
dan 



^ (V<^) • (V0) dr, 

which requires solving one forward problem to compute 
and N integrations. 



one adjoint problem to obtain 



4 Conclusion 

We presented explicit and generalizable methods of calculating the sensitivity in inversion of 
Poisson's type problems. For problems that minimize the misfit between the data and the 
model for the purpose of inversion, gradient descent methods or Newton type methods such 
as the Gauss-Newton and the Levenberg-Marquardt may be used. From an implementation 
point of view, gradient descent methods are easy to implement, but iteratively slow and 
have variable scaling issues |31|. On the other hand Newton type methods have a faster 
convergence rate and are robust to variable scalings but can be computationally expensive. 
Based on the nature of the problem and the available computing resources either methods 
may be desirable in solving an inverse problem. It is highlighted that to maintain efficiency, 
two different adjoint problems need to be solved to obtain the sensitivity information in each 
inversion scheme. This paper beside providing a step by step derivation of the sensitivities 
for the Poisson's equation with general type of boundary condition, clarifies the distinction 
and the relationship between the two forms of inversion for various types of applications. 

5 Appendix 

By definition, for the systems (22)-(23) and ([9l)-(|10[) to be adjoints we must have 



f V ■{cFV(i)i)(l)dr^ f V ■ {aV(t))(i}idT. (26) 
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Applying the identity ( 13 ) twice to /q V • {aSJ4>i)4> dr, first with g = (p and F = aV^ii, and next 
with g - (pi and F - aV(p results in 



/ V ■ {(J\/ci)i)(l)dr- / \/ ■ {cr\/ct>)(f>idr 



(27) 



We next split the surface integral on the right hand side of ( |27[ ) over regions Pq, where o. 
and r \ Fq, where a-0. On using (10) and (23) we have 



an an 



a 



and over F \ Fa that a-0 and /3 * 0, using (10) and (23) yields 
therefore 



a—(t>-a—( 
on on 



crS, d(j)i 
(3 dn 



Based on (28) and (^29| we rewrite (27) as 



/ V ■ {(J\/ci)i)(l)dr- / V • (crV0)(^i dr 



(28) 

and (j) - and 
(29) 

(30) 



f ^-^dS^ f 
Jr\ra a JVa 3 dn 



dS. 



Comparing (30) to (26) shows that for the systems ([22])-([23]) and ([9])-(10) to be adjoints 
the expression on the right hand side of (30) needs to be zero and this is not generally the 
case. Clearly (26) holds for the homogeneous case (^ = 0) and the two systems are adjoints as 
mentioned in the text. 
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